
%%
% X=10*rand(5,2);
% 
% distance=pdist(X);
% distance(distance>5)=0;% 大于20KM,断开路径可能性
% adjMatrix=squareform(distance)
% 
% judge=RDFS(adjMatrix)


%% 
clc;
clear;
close all;

load('input_1.mat');
DT=200;
location=location_NP;
upper_bound=N;
% location 输入 NP个解，分别对应 D个中继器 D*2*NP

%% Problem Definition

CostFunction=@(x) Sphere(x);        % Cost Function函数句柄，未知数是x，相当于建立了一个函数文件,类似在C语言中的函数定义。该方法在Sphere.m中定义

nVar=size(location,1);             % Number of Decision Variables问题的维度  %D维

VarSize=[ nVar 2];   % Decision Variables Matrix Size定义一个问题维度大小的矩阵

VarMin=0;         % Decision Variables Lower Bound函数的下限
VarMax=upper_bound;         % Decision Variables Upper Bound上限

%% ABC Settings

MaxIt=1000;              % Maximum Number of Iterations迭代次数上限

nPop=size(location,3);   % Population Size (Colony Size)初始雇佣蜂数量  %%%   NP个蜜源

nOnlooker=nPop;         % Number of Onlooker Bees初始观察蜂数量

L=round(0.6*nVar*nPop); % Abandonment Limit Parameter (Trial Limit)  round():四舍五入取整，表示蜜源试验次数上限，如果达到此上限，舍弃该蜜源，侦查蜂来做的这一步
% a=1;                    % Acceleration Coefficient Upper Bound蜜源变换的加速系数的最大值

%% Initialization
% Initialize Population Array

pop={};
% Initialize Best Solution Ever Found因为是求函数的最小值，所以预定义最坏适应度是无穷大
BestSol.Cost=inf;

% Create Initial Population在定义域范围内，随机初始化蜜源位置
for i=1:nPop
    pop{i}.Position=location(:,:,i);%%unifrnd(-10,10,[1 5]):表示产生均匀分布的随机数，产生一个1*5的随机数矩阵，其值在-10到10之间均匀分布
    pop{i}.Cost=CostFunction(pop{i}.Position);
    if pop{i}.Cost<=BestSol.Cost%%蜜源适应度较好的保留
       BestSol=pop{i};
    end
 end




% Abandonment Counter
C=zeros(nPop,1);%%用来记录蜜源的试验限制，达到L次数，就舍弃蜜源，用新的替代这个蜜源
% Array to Hold Best Cost Values
BestCost=zeros(MaxIt,1);%%记录迄今为止最好的蜜源的适应度

%% ABC Main Loop

for it=1:MaxIt
    
    %%%%%%%%%%%%%%%%% 雇佣蜂Employed Bees  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:nPop
      
        % Choose k randomly, not equal to i
    ith_postition=pop{i}.Position; % NP中的一个蜜源
    R(i)=size(ith_postition,1);
    dim_R=randperm(R(i));
    g=dim_R(1);
    k=dim_R(2);
        
%从 i 到 k 之间
    if g > k
         range_r=[k,g];
    else
         range_r=[g,k];
    end
    % New Bee Position 确定新一个蜜源位置，维度为R
    
   t=range_r(1):range_r(2) ;  
    temp_position=mean(ith_postition(t,:));
    temp_R=1:R(i);
    temp_R(t)=[];
    anti=temp_R;
    newbee.Position=[ith_postition(anti,:);temp_position]; 
    %检验通路是否满足要求
    distance=pdist([location_TU;newbee.Position;location_EOC]);
    distance(distance>DT)=0;% 大于20KM,断开路径可能性
    adjMatrix=squareform(distance);
    judge=RDFS(adjMatrix);
    
    if judge==0 %不满足要求，换蜜源
        continue
    end
        % Evaluation计算新蜜源的适应度值
        newbee.Cost=CostFunction(newbee.Position);
        
        
        % Comparision
        if newbee.Cost<=pop{i}.Cost
            pop{i}=newbee;
            whatt=1;
        else
            C(i)=C(i)+1;%%不好的蜜源，浏览该蜜源次数加一，浏览次数达到L次就舍弃该蜜源
            whatt=0;
        end
        
    end
  
 % 迭代结束
    % Calculate Fitness Values and Selection Probabilities
    F=zeros(nPop,1);%以下是公式去计算蜜源的选择概率
    sumCost=0;
  for i=1:nPop
       sumCost=sumCost+pop{i}.Cost; % Convert Cost to Fitness
  end
    
   for i=1:nPop
        F(i) = pop{i}.Cost/(1+sumCost); % Convert Cost to Fitness
    end
    
    P=F/sum(F);%%蜜源的选择概率
    %%%%%%%%%%%%%%%%% 雇佣蜂Employed Bees  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    
    
    
    
    
    %%%%%%%%%%%%%%%%% 观察蜂Onlooker Bees  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for m=1:nOnlooker%让每个观察蜂去干该干的事情
        i=RouletteWheelSelection(P);      
        % Select Source Site根据选择概率 进行轮盘赌选择一个蜜源
        
    ith_postition=pop{i}.Position; % NP中的一个蜜源
    R(i)=size(ith_postition,1);
    dim_R=randperm(R(i));
    
    g=dim_R(1);
    k=dim_R(2);
        
%从 i 到 k 之间
    if g > k
         range_r=[k,g];
    else
         range_r=[g,k];
    end
    % New Bee Position 确定新一个蜜源位置，维度为R
    
   t=range_r(1):range_r(2) ;  
    temp_position=mean(ith_postition(t,:));
    temp_R=1:R(i);
    temp_R(t)=[];
    anti=temp_R;
    newbee.Position=[ith_postition(anti,:);temp_position]; 
    
    %检验通路是否满足要求
    distance=pdist([location_TU;newbee.Position;location_EOC]);
    distance(distance>DT)=0;% 大于20KM,断开路径可能性
    adjMatrix=squareform(distance);
    judge=RDFS(adjMatrix);
    
    if judge==0 %不满足要求，换蜜源
        continue
    end
        % Evaluation计算新蜜源的适应度值
        newbee.Cost=CostFunction(newbee.Position);
  
        % Comparision判断蜜源留下还是舍弃
        if newbee.Cost<=pop{i}.Cost
            pop{i}=newbee;
        else
            C(i)=C(i)+1;
        end
        %%%%%%%%%%%%%%%%%%%%%%又重新更新蜜源信息%%%%%%%%%%%%
    end
    %%%%%%%%%%%%%%%%% 侦查蜂Onlooker Bees  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%% Scout Bees %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:nPop
        if C(i)>=L%%判断蜜源是否达到试验上限，达到了就丢弃蜜源
            pop{i}.Position=unifrnd(VarMin,VarMax,VarSize);%重新建立一个蜜源替代原来的蜜源
            pop{i}.Cost=CostFunction(pop{i}.Position);
            C(i)=0;
        end
    end
    
    % Update Best Solution Ever Found
    for i=1:nPop
        if pop{i}.Cost<=BestSol.Cost
            BestSol=pop{i};%%寻找最佳蜜源
        end
    end
    %%%%%%%%%%%%%%%%% Scout Bees %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Store Best Cost Ever Found
    BestCost(it)=BestSol.Cost;%%存储最佳蜜源
    
    % Display Iteration Information
    %%disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
    
end
    
%% Results
result=BestSol.Position %%打印解的结果
size(result)
%%%适应度值变换曲线图
plot(BestCost,'LineWidth',2);
semilogy(BestCost,'LineWidth',2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;





